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ABSTRACT 

A mesh deformation scheme is developed for a structured multi-block Navier-Stokes code consisting of 
two steps. The first step is a finite element solution of either user defined or automatically generated 
macro-elements. Macro-elements are hexagonal finite elements created from a subset of points from the 
full mesh. When assembled, the finite element system spans the complete flow domain. Macro-element 
moduli vary according to the distance to the nearest surface, resulting in extremely stiff elements near a 
moving surface and very pliable elements away from boundaries. Solution of the finite element system for 
the imposed boundary deflections generally produces smoothly varying nodal deflections. The manner in 
which distance to the nearest surface has been found to critically influence the quality of the element 
deformation. The second step is a transfinite interpolation which distributes the macro-element nodal 
deflections to the remaining fluid mesh points. The scheme is demonstrated for several two-dimensional 
applications. 


INTRODUCTION 

The simulation of complex flow fields such as vortex, separated and transonic flows requires nonlinear 
computational fluid dynamics (CFD). Modern multi-discipline CFD often also incorporates moving boundaries, as 
would be required in applications such as design optimization, aeroelasticity or forced boundary motion. These in 
turn require an algorithm for the deformation of the mesh. Considerable effort has been devoted to the development 
of robust and efficient general purpose techniques for mesh deformation. And as techniques improve, the difficulty 
and complexity of the applications has advanced. The aeroelastic simulation of flutter and buffet response of multi- 
engine transports and the large amplitude deflection of highly flexible large aspect ratio and novel aircraft 
configurations are at the forefront of current challenges. New challenges will yet evolve with greater complexity of 
motion. For these reasons continued development of better methods of mesh movement is worthwhile. 

Mesh deformation schemes can be separated into algebraic methods, elasticity based methods, or a combination of 
the two. Algebraic methods would include grid shearing, transfinite interpolation (TFI) with various blending 
functions, exponential decay of surface deformation, or a combination of shearing and rotation. References 1 and 2 
offer algebraic methods in which surface movement is transmitted with an exponential decay into the mesh interior 
via selected control points. Intervening mesh points are then moved by a TFI step using arc length blending. 
Approaches based by analogy on elasticity would include the spring analogy with tensile (also called lineal) and/or 
torsion spring stiffness, and those methods based directly on the equations of continuum mechanics. The spring 
analogy scheme for mesh movement, first introduced by Nakahashi and Deiwert (ref. 3) and later adapted for 
aeroelastic applications by Batina (ref. 4), has been the mainstay technique for mesh deformation owing to its 
simplicity and efficiency. It is well known, however, that the standard spring analogy can result in mesh 
irregularities especially for large deformations. Various attempts to enhance the robustness of the spring analogy 
without substantially increasing the computational effort, notably the stiffening of springs near moving boundaries 
and semi-torsion springs, have resulted in some success, (ref. 5) The introduction of nonlinear torsion stiffness has 
been shown in reference 6 to significantly enhance the robustness of the scheme. This also increases somewhat the 
computational overhead required. Algebraic and elasticity based approaches have also been combined to form 
hybrid approaches. Hybrid approaches have been constructed with surface translation and rotation transmitted into 
the near field region decaying away from surfaces to a mesh motion computed using the spring analogy, (refs. 7-8) 

Other elasticity based schemes have started from the equilibrium equations of continuum mechanics. Reference 9 
constructs a solution for mesh deformation in which layers of near field mesh points move rigidly while all other 
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mesh points are moved by the finite element method (FEM). Smaller elements are made stiffer relative to larger 
ones by dropping the Jacobian of the transformation from the element to physical domain. In reference 10 the 
mesh motion obeys the equations of isotropic linear elasticity. A spatially varying Poisson’s ratio is employed in the 
attempt to maintain mesh integrity for highly skewed cells. Material behavior is based on aspect ratio and size of 
elements. This approach has been shown to perform better than the original tensile spring analogy. Reference 1 1 
uses a boundary element solution of the equations of elasticity with nonlinear stiffness based on local grid volume to 
deform the mesh. 

Each of the foregoing approaches has solved the elasticity equations using all or part of the mesh used to solve the 
fluid equations. This has the advantage of simplicity in that the fluid mesh becomes the de facto finite element 
mesh. On the other hand, it does require the solution of a large finite element system. Solution of the mesh 
deformation can easily require the same or more solution time as the flow field. It also has the disadvantage that the 
fluid mesh clustering has been designed with the flow field in mind rather than with a consideration of the mesh 
deformation to be resolved. 

In this paper an approach is presented that performs a finite element solution to obtain mesh deformation in a 
structured multi-block finite volume CFD code. Since the CFD solver is a finite volume based procedure, the finite 
element method is the most appropriate technique for deforming the mesh. The present mesh deformation separates 
the finite element solution of a selected number of points from the entire flow field mesh. The deformation of 
selected points is transmitted to the entire mesh with a TFI step. There are several important features to the present 
approach. The use of a subset of control points, also called macro-element node points in this paper, reduces the 
computational overhead required for the finite element solution. Another advantage is that the use of control points 
allows the user to reduce the disparity of element sizes and thus the numerical stiffness that will be associated with a 
high Reynolds number grid. Finally, although this approach is most easily applied to a structured code, it should be 
possible to extend the ideas presented here to an unstructured mesh. The next section presents the approach in detail, 
while a subsequent section applies the approach to several two-dimensional deforming meshes. 

APPROACH 

Finite element based mesh deformation approaches have typically used the entire fluid flow mesh in defining the 
finite elements. In the present approach a finite element solution is performed as the first of two steps in the process 
of deforming the fluid mesh. The finite element model is constructed with a subset of mesh points called macro- 
element node points which are element vertices. Figure 1 presents the orientation, coordinate system and vertices, 
or node points, for a typical macro-element in the computational and physical domains. In a structured grid code the 
step of identifying element vertices is relatively straight forward. 

Figure 2 shows several ways that nodes can be specified. One way is to specify skip values for each index direction 
and block. Figure 2a shows node points spaced with constant mesh skip values in the three index directions. 

Another way to specify nodes is to arbitrarily select block and index locations, as shown in Figure 2b. In that figure, 
physical sizing of macro-elements has been made more uniform by varying the number of mesh points skipped from 
one macro-element to the next. The flexibility of specifying arbitrarily located node points allows the user to size 
the macro-elements as needed. Arbitrary placement further allows the user to define the smallest number of macro- 
elements necessary to adequately resolve the anticipated surface deformation, thus reducing the computational effort 
required. The reduced computation and memory requirement also allow the identification and storage of element 
connectivities at the start without significantly increasing the memory requirement of the overall code. The only 
constraint is that all block boundaries must be included as macro-element faces. Block boundaries as macro-element 
faces results in coincident points at blocking interfaces and adds a minor complication to the assembly process. 
Finally, the deflections obtained from the finite element solution using these macro-elements provides the 
displacement of the subset of the element vertices. The second step in the deformation of the flow field mesh is the 
interpolation of deflections between macro-element vertices using transfinite interpolation with an arc length 
blending function. This step will be taken up later. The finite element solution step will be discussed in greater 
detail first. 

The finite element solution follows the standard formulation beginning with the definition of material properties. 

The constitutive relationships between stresses and strains are given by Hooke’s law for a macro-element m. 


2 



where 


o„ 


C m £ m 


°xx 


£ xx~\ 


1 

3 

0 

0 

0 

0 

0 “ 

°yy 


e yy 


0 

E , n 

0 

0 

0 

0 

O 


£ „ 


0 

0 

E m 

0 

0 

0 

zz 

, £„ = 

ZZ 

, c = 


m 




°xy 

5 m 

£ xy 

5 m 

0 

0 

0 

G m 

0 

0 

°yz 


£ yz 


0 

0 

0 

0 

G m 

0 

1 

to 

1 

m 

Jxz\ 

m 

0 

0 

0 

0 

0 

1 

5 

o 


Since induced lateral strains do not materially contribute to the robustness of the mesh movement, these are set to 
zero (i.e. Poisson’s ratio, v= 0). The macro-element moduli are constant throughout a given element, and are 
modeled by a fictitious isotropic inhomogeneous material. In the context of the spring analogy scheme, reference 5 
shows that there is an improved robustness for large deflections by the simple expedient of increasing the spring 
stiffness near the moving boundary. Along that vein, the present finite element based approach varies material 
properties of the macro-element (referred here after simply as “material properties’’) according to the distance of 
elements from the nearest mesh boundary. A function exhibiting an exponential rise with decreasing distance from 
a boundary is the basis for the magnitude of material properties. For the element m , the Young’s and shear moduli 
are obtained by 
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where the subscript ( ) mc denotes the center of element m. The reference Young’s modulus E 0 is arbitrarily 
specified. The parameter /? controls the rate of decay of material properties away from surfaces. Typical values are 
0(1 ). The subscript ( ) lm denotes the nearest surface face center or nearest surface node point which ever results in 
a smaller value of A r . Figure 3 illustrates this distance for several representative elements. Note that for element 

m 1 the nearest surface point to the center of the element is the point sj located at the center of the surface face. For 
the element m 2 located just beyond the trailing edge the nearest surface point is the node point s 2 at the trailing edge. 
The definition of distance is crucial to ensure that element moduli vary smoothly and have appropriate magnitude 
near moving surfaces. Although not entirely general, the current approach has resulted in fairly smooth variations in 
mesh deformation. This feature will be explored in the examples to follow. The choice of the function f m in equation 

(1) for the moduli decay is determined by the fact that it has the following properties: as A r m — > 0, f — * 

and as A r m — > oo, These properties produce a material stiffness approaching infinity as elements are 

compressed into and become coincident with a moving surface. Although material property variation as a function 
of mesh deformation results in an inherently nonlinear problem, the material properties for each element are 
computed once per time step. Computation of properties at each time step ensures that as elements compress toward 
or stretch away from moving surfaces local material stiffness increases or decreases accordingly. 
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Figure 1 indicates the local coordinate system £,^>77 of a hexagonal macro-element. The stiffness matrix for element 
m is obtained from 


1 1 1 

k = \b t C B cIV = J- x \\\b t C B d£d£dn 

m J m m m rn III m m m ^ ^ / 

v m 0 0 0 

where B,„ relates the six stresses to the 24 nodal displacements of element m, i.e. 


= B.U„ 



J is the Jacobian of the coordinate transformation and N m is the subset of node points that form the vertices of 
element m. The elements are assembled into a global stiffness matrix and the linearized global finite element 
equation 


KU =R 


is solved. The solution vector U contains the N nodal displacements 
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where U a is the vector of off-surface nodal displacements, U s is the vector of surface nodal displacements. Surface 
nodal displacements are prescribed by dividing the stiffness matrix K and the right hand side R into surface and off- 
surface partitions. 


K = 


K° K° 
0 1 


R = 



where (U s ) p is the prescribed surface displacement. The entire asymmetric system is solved by the bi-conjugate 
gradient minimum residual (BCGMR) method with diagonal element preconditioning. 

Once the macro-element node displacements and updated grid locations are obtained, the second step is performed, 
namely the computation via transfinite interpolation of the displacements of interior mesh points. This is a natural 
way to transmit nodal displacements to all interior mesh points since the strains used in derivation of the stiffness 
matrix are also derived based on a transfinite interpolated representation of deflections within the element. The 
result is a continuity of displacements at all macro-element boundaries, although continuity of the derivative of the 
displacement is not ensured. This feature will be explored in the next examples. 
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EXAMPLES 


Pitching Airfoil 

The first example consists of a pitching NACA 0012 airfoil. Figures 4 and 5 show a viscous grid about the airfoil 
and the corresponding user-defined macro-elements. In the present example, since the airfoil motion is pitching, a 
coarse and fairly uniform macro-element distribution has been specified. Figures 6-7 present the final macro- 
element orientations and mesh after the airfoil is pitched to 50 degrees. It can be seen in Figure 6 that the effect of 
the displacement of the airfoil trailing edge is transmitted a considerable distance down stream from the airfoil. This 
effect manifests the influence of the element stiffness decay away from the airfoil. For this mesh the decay 
parameter ft is set at 6. The mesh motion for all the cases studied appears to be fairly insensitive to parameter values 
at least in the range 1 < < 6 . The nearly rigid rotation and translation of macro-elements near the airfoil 

surface, as seen in Figure 7, is evidence of the large element stiffness near the airfoil surface. This feature results in 
the preservation of the mesh orientation near the airfoil. The coarse macro-element definition of this example 
produces the appearance of facets in the final fluid mesh (apparent in Figure 7b) due to slope discontinuities in the 
mesh at junctions between elements. A finer macro-element distribution near the airfoil surface would reduce the 
prominence of this feature. 

The robustness and sensitivity of the scheme to orientation and spacing of macro-elements is an important 
consideration. To test sensitivity to element spacing, an alternate definition of the macro-elements for the pitching 
airfoil are shown in Figure 8. The macro-elements are defined identically to the previous case, except that the first 
stream wise line of nodes off the surface are spaced more closely to the airfoil. As seen in Figure 9, the elements 
near the airfoil surface retain shape well after rotation of the airfoil to 50 degrees. There is a slight flattening of the 
leading edge element which will result in some reduction in grid quality near the leading edge. Another feature can 
be seen in Figure 9b, by observing the element side proceeding vertically upward from the airfoil trailing edge. That 
line does not bisect the angle created by the airfoil upper surface and the node line proceeding down stream from the 
trailing edge, but rather remains nearly normal to the airfoil upper trailing edge surface. While mesh points over the 
upper surface retain their original orientation relative to the airfoil surface, an undesirable consequence is an 
excessive compression of mesh points beyond the trailing edge. More appropriate mesh motion would result from 
an even compression of all elements ahead of and behind the trailing edge. 

The reason for deformation of the macro-elements in Figure 9b near the trailing edge is evident. Because the 
elements running along the airfoil surface are extremely stiff, they deform very little. The elements proceeding 
down stream from the trailing edge have centers far from the airfoil surface, and thus are quite flexible. This 
disparity in element moduli results in the uneven element distortion. An intuitive approach to this disparity would 
be to add vertical node lines beyond the trailing edge to create smaller macro-elements in that region that will have 
centers much closer to the airfoil surface, and thus reduce flexibility near the trailing edge. Such a distribution of 
elements results in the solution shown in Figure 10. By comparing the solutions in Figures 10 and 9b it is clear that 
the element just behind the trailing edge is not distorted quite as much in Figure 10 as it is in Figure 9b. However, 
the compression of the elements ahead of and behind the trailing edge is still not as smooth as would be desired. 

That some distortion remains, and will most likely always be present to some degree, can be explained by referring 
again to the distances shown in Figure 3. Because the distance from the face center on the airfoil surface to the 
center of the element m I is smaller than the distance from the trailing edge to the center of the element m 2 , element 
m 2 will always have less stiffness than element m 2 . This fact implies that the present method of defining distances 
for computation of element moduli will always result in some distortion of elements at a surface discontinuity such 
as a trailing edge. It also illustrates how critical the computation of element moduli is to the smooth deformation of 
a mesh. 


Multi-element airfoil with flap rotation 

The second example illustrates use of the scheme for a complex multi-block viscous grid about a three-element 
airfoil. The three-element airfoil is composed of a leading edge slat with a leading edge down 30 degree deflection, 
the main airfoil section and a trailing edge flap, initially deflected 30 degrees trailing edge down. The flow field 
mesh is composed of three C-grid blocks enclosing each of the airfoil sections and a fourth C-grid enclosing all the 
sections and grids. The fluid mesh has a total of 271,000 grid points. To set up the required mesh deformation, 664 
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points distributed over the four flow field blocks, are identified as the nodes of the macro-elements. A time accurate 
simulation is then performed with the trailing edge flap rotating downward about the point X = 0.72c, y = 0 
(which happens to be forward of the flap leading edge) to an additional 30 degrees in about 100 time steps. This 
additional rotation results in a rotation and downward translation of the flap. 

The results of this simulation are found in Figures 1 1-18. The initial mesh about all three elements (with the flap at 
a 30 degree deflection) and that for a 60 degree flap deflection are shown in Figure 11. Figure 12 presents the initial 
30 degree and 60 degree macro-element orientations. With the exception of the flap leading edge, elements near the 
flap rotate and displace almost rigidly with the flap. As seen in Figure 13, the stiffness of macro-elements near the 
flap maintains the relative mesh orientation near the flap surface and overall mesh spacing over a large region of the 
wake downstream of the main wing-element. While the macro-elements apart from other constraints are free to 
move nearly rigidly, elements within the gap between the main section trailing edge and flap leading edge are forced 
to strain. Figures 12 and 13 reveal a smooth continuous stretching of elements within that region. The continuity of 
strains between macro-elements, enforced by the finite element solution, results in an equally smooth fluid mesh 
stretching. 

The efficiency of the solution can be assessed by considering the convergence rate of the BCGMR scheme and 
computational overhead required for this example. The FEM in this example has 664 node points and a total of 1992 
degrees of freedom. The BCGMR algorithm requires for this example typically about 400 iterations to converge the 

residual to 1 X 10 9 . For a flow solution taking five sub-iterations per time step, the movement of the mesh and re- 
computation of metrics requires an additional 13 percent of CPU time. 

When using the spring analogy scheme, the nonlinearity of the spring stiffness can at times be a cause of mesh 
irregularities. While the stiffening of springs as mesh points approach one another precludes grid crossing in the 
spring analogy. On the other hand, the weakening of spring stiffness as they move apart can be the cause of 
anomalous isolated gaps in the mesh (ref. 10). The results of Figures 12 and 14 show a smooth variation in 
displacement of nodes in the flap cove area, with no tendency for anomalous gaps or unusually large localized 
stretching of one element compared to the next. Another weakness of the spring analogy scheme is the lack of 
robustness of the mesh movement at discontinuities such as a trailing edge. In the present example, because of the 
stretching of the elements just behind and below the main wing-element trailing edge as the flap moves downward, 
large gradients in grid deformation occur. In such a situation, the spring analogy scheme can potentially pull the 
tightly spaced upper surface grid lines down over the trailing edge, resulting in a grid collapse (ref. 10). The present 
results show that the macro-elements above the trailing edge of the main airfoil section retain nearly their original 
shape and size, thus maintaining mesh integrity over the top trailing edge of the airfoil. 

Figure 15 shows the return of the flap element beyond its original position to a deflection less than 30 degrees. The 
flap leading edge is displaced upward into the main wing-element trailing edge. This results in the collapse of 
macro-elements that is seen in Figure 16 in a localized region near the trailing edge of the main wing-element. 

Figures 17 and 18 show the same flap motion for macro-elements that have been automatically generated by a 
meshing algorithm. The algorithm creates the fewest number of macro-elements possible under the constraint of 
constant i,j and k index skip values. In the present example the constraint of constant index skip values forces the 
creation of many more macro-elements than spaced elements. Although all macro-elements have positive volumes 
and no mesh irregularities occur, the appearance of overlapping grid lines below the leading edge of the fully 
deflected flap hints at the potential approach of grid collapse. A solution to this potential collapse is to define an 
additional stream wise line of nodes below the flap, as was done with the macro-elements of Figures 8-13, forcing a 
more regular spacing of elements. 


CONCLUDING REMARKS 

A mesh deformation scheme has been developed using macro-elements composed of node points selected from the 
fluid mesh. A finite element solution using these macro-elements updates the node point mesh displacements. A 
fictitious isotropic, inhomogeneous material is used in which stiffness moduli computed with an exponential 
function based on distance from the nearest mesh boundary. As the distance of a macro-element center to a moving 
surface tends to zero, the moduli grow rapidly. This growth results in very stiff elements near a moving surface and 
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very pliable elements away from boundaries. A second step in the scheme interpolates these displacements to the 
remaining mesh using transfinite interpolation. This approach has been shown to have several advantages as well as 
a disadvantage compared to other approaches. Its disadvantage is that additional effort is required in the initial 
definition of the macro-elements. On the other hand, the choice of node points for macro-element vertices allows the 
judicious placement of elements as needed to resolve the anticipated displacement. The freedom of node placement 
contributes to a potentially dramatic reduction in computation time compared to a finite element solution over the 
entire mesh and allows asymmetries in mesh sizes to be reduced in the macro-element definition, thus enhancing 
robustness. 

The use of highly stiff elements near moving surfaces has been shown to result in nearly rigid translation and 
rotation of those elements with the surface motion. This characteristic results in the maintenance of high mesh 
quality in areas of large deformation. The scheme does, however, tend to produce slope discontinuities in the mesh 
at macro-element boundaries. This feature can be alleviated by adding more intervening macro-elements to 
distribute the deformation more smoothly. The quality of element deformation has been found to be quite sensitive 
to the manner in which distances are computed. A distance function has been demonstrated that produces a 
relatively smooth mesh deformation. 

The current work has also shown the results of an automatic meshing technique that reduces user effort in defining 
macro-elements. The algorithm has been tested for a complex 2D multi-block mesh and found to produce 
acceptable mesh motion. Future work will require that this scheme be validated for three-dimensional cases and 
implemented on a parallel computing environment. The method is naturally suited for the wide variety of topologies 
seen in three-dimensional finite volume grids for complex geometries. Mapping of macro-elements to multiple 
processors in the most efficient way will be required, as well as the parallelization of the assembly and solution of 
the linear system. Work will also be required to assess the robustness and sensitivity of the scheme to a wide range 
of geometries and element orientations. 
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a) Leading edge macro-elements 


b) Trailing edge macro-elements 


Figure 9. Leading and trailing edge detail (alternate macro-element definition), 
airfoil pitched to 50 degrees. 



Figure 10. Trailing edge detail (alternate macro-element definition), 
airfoil pitched to 50 degrees. 





a) Initial mesh, flap 30 degrees 


b) Final mesh, flap 60 degrees 


Figure 1 1 . Three element airfoil with deflecting flap. 


a) Initial macro-element orientations, flap 30 degrees 


b) Final macro-element orientations, flap 60 degrees 


Figure 12. Macro-elements for the three element airfoil. 


a) Initial mesh, flap 30 degrees 


b) Final mesh, flap 60 degrees 


Figure 13. Trailing edge and flap mesh detail. 
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Figure 17. Automatically generated macro-elements, constant skip, flap 30 degrees. 



Figure 18. Automatically generated macro-elements, constant skip, flap 60 degrees. 
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